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Abstract 

We consider the drift and diffusion properties of periodically 
driven renewal processes. These processes are defined by a period- 
ically time dependent waiting time distribution, which governs the 
interval between subsequent events. We show that the growth of the 
cumulants of the number of events is asymptotically periodic and 
develop a theory which relates these periodic growth coefficients to 
the waiting time distribution defining the periodic renewal process. 
The first two coefficients, which are the mean frequency and effective 
diffusion coefficient of the number of events are considered in greater 
detail. They may be used to quantify stochastic synchronization. 



1 Introduction 

Many dynamical processes in physics, biology and chemistry, although be- 
ing inherently continuous, can be reduced to a series of discrete events 
without losing much information. Often this reduced description in terms 
of discrete events has the additional property, that time intervals between 
subsequent events are statistically independent. Such processes are called 
renewal processes IJ. They are fully described by a waiting time distribu- 
tion w{t), which governs the statistical properties of the time intervals r 
between two subsequent events. In the simplest case the probability per 
unit time, i.e. the rate, for an event to occur is independent on the time 
elapsed since the last event. This subclass of renewal processes is called 
Markovian, because the event number of such a process as a function of 
time is a Markovian stochastic process. If driven externally, the rate can 
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be temporally modulated becoming thus a function of the absolute time 
t, but still remains independent on the waiting time t elapsed since the 
previous event. 

In the general case, the probability per unit time that an event will 
happen is not independent of the time already elapsed since the last event. 
Then a non exponentially distributed waiting time between subsequent 
events is observed which distinguishes these general renewal processes from 
discrete Markovian processes. 

The concepts of renewal processes have been successfully appHed to 
random sequences of spiking events in the dynamics of neurons IH E] , 
to random walks in a tilted periodic potential |2I and to the failure times 
of biological and technical machineries [J]. In these models the processes 
which generate the single events consist of a series of several Markovian 
steps or possess a priori compHcated waiting time densities. For example, 
the generation of a new spike requires excitation of the voltage variable 
over a threshold value which is followed by the spiking and refractory time 
where neurons are unable to produce new spikes. Also the ionic trans- 
port through pores involves several steps as diffusional motion followed by 
escapes over potential barriers. Molecular motors walk along microtubuH 
where sequences of different configurational changes of the proteins realize 
one forward step. 

There is a long history of studies of stationary renewal processes QP, i.e. 
renewal processes having waiting time densities w{t), which do not depend 
on the absolute time t. Also stationary situations with multiple events and 
different stationary waiting time statistics were studied. Quantities of in- 
terest like the mean event number in a certain period of time or the event 
number diffusion are known in terms of the waiting time distribution, be- 
cause a fully developed theory exists for these time homogeneous processes 

mm 

In connection with stochastic resonance (SR) the interest for non sta- 
tionary but periodic processes has grown [HliniEl- In particular, Marko- 
vian periodically driven models have been investigated ^I]. The two state 
Markovian theory with non stationary rates describes successfully the dy- 
namic behavior in bistable situations as shown in many different studies in 
various fields of science. However, other systems Hke excitable dynamics as 
used for example to describe the spiking mechanism in neuronal systems 
can no longer be approximated by a simple discrete Markovian description. 
If such systems are periodically driven, as is for example the case in neu- 
rons which respond to periodically varying inputs, more general concepts 
are needed to describe their behavior. In this paper we study the general 
situation of non stationary but periodic renewal processes with arbitrary 
waiting time distributions w{t, t). The periodicity is reflected by a periodic 
dependence of the waiting time distribution on the absolute time t of the 
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previous event. 

There exist different possibilities to quantify the periodicity of a peri- 
odic stochastic process, or, if one considers this periodicity as induced by 
an external periodic signal, the quality of the response of the system to this 
periodic signal. On the one hand, spectral based measures like the spec- 
tral power amplification and the signal to noise ratio have been frequently 
employed (see jlOj and references therein). These spectral based response 
measures where also considered in the context of periodically driven renewal 
processes in El and for a special discrete state model for excitable dy- 
namics in pHl- To this end the sequence of events has to be somehow 
mapped onto a stochastic process, e.g. by considering a sequence of delta 
peaks located at the event times or by assigning alternatingly after each 
event different values to the process. Another possibility to characterize 
periodicity of the process is to consider the drift and diffusion properties 
of the number of events. The evolution of the number of events in time 
can be characterized by two quantities, the mean frequency of events and 
the effective diffusion coefficient, which describes, how the variance of the 
number of events grows in time and thus characterizes the regularity of 
the process. The lower this effective diffusion coefficient for a fixed mean 
frequency, the more regular, i.e. periodic, is the system's dynamics. Both 
quantities together may therefore serve as a measure of periodicity of the 
process and thus, if this periodicity stems form the infiuence of a periodic 
input, as a measure of stochastic synchronization between the driving signal 
and the system dynamics [HI [121 ITfil [131 [21] . 

In this paper we present a method to calculate the mean frequency and 
the effective diffusion coefficient of a periodically driven renewal process in 
terms of its periodically time dependent waiting time distribution w{T,t). 
After having derived the general concepts and results in sections 12 to 
we consider two simple situations, namely general but undriven renewal 
processes and periodically modulated rate process, for which an expHcit 
evaluation of the general results is possible and agrees with the known 
results. Finally in section and |H1 we numerically evaluate our theory for a 
toy model where the intervals between subsequent events are governed by 
a a fixed but periodically varying waiting time followed by a rate process 
with constant rate. The results agree with simulations of the underlying 
periodic renewal process. Many of the calculations are summarized in a 
series of appendices. 

2 Periodically driven renewal processes 

A periodically driven renewal process is a sequence of events occurring at 
times . . . , ii, ii+i, ti+2, ■ • ■■ The intervals between two subsequent events 
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are governed by the waiting time distribution w{T,t). The first argument 
r represents the waiting time whereas the second argument t denotes the 
absolute time at which the last event happened. Therefore, w{t, ti)dT is the 
probability that the event i+l happens in the time interval (ti+r, ii+r+dr] 
if event i have been happened at time ti . The dependence of the waiting 
time distribution u'(t, t) on time t is due to the periodic driving and thus 
periodic with the period T = 2t: /Vt of the signal. Although, in contrast to 
an ordinary renewal process, the intervals between subsequent events are 
now correlated, it is still only the time of the previous event which governs 
the statistics of the following event, which justifies to stick to the term 
renewal. Normalization holds at arbitrary time t 

dTw{T, t) = 1 



An) , ^u,,+At ^to, ^ a i ) 





A lot of information about these processes is contained in the random 
number of events -/Vto,t which take place in the interval (io, t]- For example 
stochastic synchronization to the periodic driving can be characterized as 
an integer relation between driving frequency and the frequency of events 
and at the same time, a decrease in the effective diffusion coefficient of the 
events, i.e. a more regular (periodic) behavior. 

We will evaluate the mean frequency and effective diffusion coefficient 
based on the periodically time dependent waiting time distribution w{t, t), 
which defines the periodically driven renewal process. To this end we con- 
sider more generally the nth cumulants K^^ ). of the event number Ntg^t and 
define their increment per time 

'(") _ T^in) 

■ta,t+At ^^tp.t 

The first of these coefficients is the mean frequency while the second is the 
effective diffusion coefficient. We show that asymptotically 

«(")(i):= lim 4:.) (2) 

become periodic functions of time with the period of the external driving. 

To this end we introduce following [2] the generating functional 
Lta.tW\ of the considered driven renewal process as 

UoAv] = (T[{l + v{U))\ 

where ti are the times of the events in the interval (to,t]. Then the nth 
moment M^^ j. of event number Nt^^t is given by 

m/" := (iV" ,) = Je" - 11 . (3) 
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The generating functional Lt^^t can be expressed in terms of the distribution 
functions fs{ti,. . . ,ts), which govern the probabiHty 

dP = fsiti,...,t,)dh...dts 

to find one event in each of the intervals {ti,ti + dti), i = 1, . . . , s regardless 
of how many events are outside these intervals, as |2] 



— 1 to -^^0 



dTs fs{Tl,T2, ...,Ts) v{ti) . . . v{Ts) 



(4) 



The generating functional can also be expressed in terms of the correlation 
functions gs{ti, . . . ,ts) as 



ito.tH = exp 



to 'J to 



dTi 



dTs gs{Tl,T2, ■ ■ ■ ,Ts)v{ti) . . .v{ts) . (5) 



Eq. lOl together with eq. Q define the correlation functions in terms of 
the distribution functions. 

According to eq. Q the moments eq. (pj can be expressed as 



(n) 
to,t 



^exp[^G.(to,t)(e"-ir 



«=o 



(6) 



where 



Gsito,t) := / dTi / dT2 . 



dTsgs{Tl,T2,...,Ts). 



From formula ^ we can evaluate the corresponding cumulants -R^t"^^ as 
(see appendix El 



ti=0 



(7) 



As the considered renewal processes are periodic in time with period T, 
the distribution functions and therefore also the correlation functions are 
Hkewise periodic in time, 



5s(ti, . . . ,Ts) = 5s(ti +T,...,Ts+T). 



(8) 
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Then the time derivative of the function Gs{ta,t) yields 



Gs{to,t) = / dT2... dTsgs{t,T2, . . . ,Ts) 

Jto Jto 



dt 



1 



1 



dT2 . . . / dTsgs{t,T2, ■■■,Ts) 



- 1)! Jo 

which can be expressed in the asymptotic Umit as 



dT2.. . / dTsgs{t, t - T2, . . . ,t - Ts) 



d 1 /"GO f-OC 

lim —Gs{to,t) = — / dT2... dTsgs{t, t - T2, . . . ,t - Tg) . (9) 



To ensure that this Hmit exists, we additionally suppose that gs{Ti, . . . , Ts) 
decreases sufficiently fast to zero for any pair of time difference In—Tj] — s- 
oo. In j2] this property is called cluster property, while for stationary sys- 
tems this property is called ergodicity. 

According to eq. Q the asymptotic time derivative eq. © is a periodic 
function in t and thus (cf. eq. 10) the coefficients 

lim (10) 

to^-oo dt 

are periodic in time, as well. 

The time dependent waiting time distribution can be expressed in terms 
of the distribution functions as pj 

w(T,t) =/l(t + r)+ V^-^ / dti...dt, fs+l{t + T,ti,...,t,). 

However, we are faced with the inverse problem, namely to express the 
distribution functions fs or Hkewise the correlation functions in terms 
of the waiting time distribution w(r, t) in order to finally evaluate the 
coefficients K^^\t) according to eqs. ifTTHl and iQ. Such a relation is not 
known to the authors and even if it exists, the explicit evaluation of the 
infinite sum in eq. 10 will be challenging. Thus to obtain the n'^^^t) in 
terms of the time dependent waiting time distribution w{T,t) we have to 
adopt a different approach. 



3 The microscopic master equation 

Let us consider the probabilities Pk{t) to have had k events up to time t. 
Furthermore let jk (t) be the probability fiux from state k to state fc + 1 , 
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i.e. the probability per time that the k + 1^* event happens at time t. This 
probabiHty obeys the continuity equation 

j^Pkit) - Jk-i{t) - Jkit). (11) 

If we further assume as initial condition that event 1 happened at time 
to, i.e. jo{t) — S{t — to), the relation between the probability fluxes of the 
renewal process can be expressed by the "microscopic" dynamics as (for the 
undriven case see e.g. j2(^ l 

jfe(i) = / dt'jk-i{t')w{t-t\t'), k>l. (12) 

J to 

Using this relation one readily obtains from the continuity equation ltTT|l 

Pk{t) = f dt'jk-iit')zit~t',t'), k>l (13) 

J to 

where z{T,t) = 1 — dr'tulr' ,t) is the probability to wait longer than r 
until the next event, if the last event happened at t. In case of a Markovian 
renewal process with time dependent rate ^{t) the probability flux jk is 
related to the probability pt by jk{t) = "f{t)pk{t). Thus in this case the 
dynamics can be completely expressed in terms of the probabilities pk- In 
the general case however we need a formulation in terms of pk and jk as 
expressed in eqs. Ijl3ll and II12II . 

The moments of the number of events Nto,t in the interval {to,t] can be 
expressed in terms of the pk as 

oo 

Mi:l^Y.^>k{t). (14) 

fe=0 

Note that the moments differ from the moments deflned by eq. ^ 

since the corresponding event number iVj^ t is conditioned on having had 
an event at time to in contrast to -/Vt^ t. However the asymptotic behavior 
of both families of moments agrees. 

In principle one can calculate the cumulants from the moments eq. Ijl4|l 
to obtain eventually the coefficients ^("^(t) according to eq. ifTTijl . However, 
in practice this is not feasible in general, as one has to calculate an infinite 
sum over the Pk{t) where each pk-, according to eqs. and Ijl2|l . is 

a /c-fold integral involving the waiting time distributions w{T,t) and the 
corresponding survival probabilities z(r, t). 



7 



4 Embedding the discrete dynamics into a con- 
tinuous one 



To find a simpler relation between the periodic coefficients K'^^\t) and 
the time dependent waiting time wlr^t), which governs the microscopic 
dynamics, we construct a continuous embedding in the asymptotic limit 
to — > — oo. Consider the envelope V{x,t) of the discrete probabilities Pk{i) 
as a probability distribution on a continuous state space. Thus, respecting 
the normalization we adopt the relation (cf. Fig. 

Pk{t)^ / dxV{x,t). (15) 

Assuming that the variation of V{x,t) within an interval (fc — i, fc + i) 




-• 1 • 1 • 

k fe+i k+1 fc+l+j k+2 



Figure 1: Schematic view of the discrete event dynamics Pkit) and the 
continuous description in terms of x. 

becomes arbitrary small for large t the difference between the moments 
and therefore also the cumulants of both the discrete and the continuous 
process tends to zero in this asymptotic limit. Accordingly, also the change 
in time of the cumulants of the continuous process have to agree with the 
corresponding quantities n^^^t) of the discrete process. 

One possible choice to formulate a dynamics of the continuous envelope 
V{x,t) is a Kramers-Moyal equation. This choice renders the underlying 
stochastic process x{t) Markovian. To ensure the equality of the increments 
of the cumulants, the coefficients of the Kramers-Moyal expansion have to 
coincide with the coefScients k^"^ (t) of the renewal process (cf. appendix 
O d), i.e. 

r\ OO r '\\tl f)n 

a^^(-'^)-E^'^^"H^)^nx,t) (16) 
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The probability current jkit) in the discrete system is correspondingly 
related to the probability current of the continuous envelope description 
J{x,t). According to the relation between the discrete and continuous 
probability eq. ltT5|l . the discrete probability current jk{t) from fc to fc + 1 
is equal to the continuous probability current J^{x,t) at x = k + ^, 

jk{t)=J{k + ^,t). (17) 

The continuous probability current J{x,t) is related to the probability 
distribution V{x,t) by the continuity equation 

l-Vi^.,t) ^ ^§-JM. (18) 
and therefore according to eq. ifTfiji 

n— 1 

Thus from eq. ifTTjl we deduce 



Finally we want to mention that the more general embedding of the 
discrete process 

rk+l 

Pk{t) ^ / dxr{x,t) and jk{t) ^ J{k + l,t), Z arbitrary. (21) 

Jk-l+l 

leads to the same results as the embedding chosen in eqs. II15II and Ill7p . 
which corresponds to choosing ? = i in eq. 1I2TI1 . 



5 The asymptotic drift and the diffusion coef- 
ficient 

Having fixed the relation between the probabilities and probability fluxes 
of the discrete renewal process and the continuous embedding, it is now 
possible to relate the coefficients K'^^\t) appearing in the continuous de- 
scription Hl(i|l to the waiting time distribution w(r, t) of the renewal process, 
involved in the microscopic dynamics Ijl2|l and Ijl3|l . As we are considering 
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the asymptotic behavior we have to pass to the asymptotic Hmit in eqs. 
Ijl3|l and Ijl2|l by shifting the initial time to —oo. This results in 



Pk{t)= dTJk-l{t-T)z{T,t-T) (22) 

JO 

and 

/•oo 

jk{t)^ dTjk-l{t-T)w{T,t-T). (23) 

Inserting eqs. lfT5|l and ll2?Hl into the above eq. ll22ll we end up with 

/•oo 

dAxP{x~Ax,t) = - / dTz{T,t-T) (24) 

Jo 

y iziL«(")(t _ r)^^V{x -]-,t-T) 



dx" 

n—1 

with x = k. The probability 'P{x — Ax,t — r) can be expressed in terms 
of the probability V{x,t) and its derivatives ^^V{x,t) by performing a 
Taylor expansion of 'P{x — Ax, t — t) around x, t and converting the time 
derivatives to derivatives with respect to the state x using the Kramers- 
Moyal equation Ijl6|l . This results in (cf. appendix EJ 



V{x-Ax,t^T) = V{x,t) + c[^\t-T,Ax)—r{x,t) + 0{2) 

ox 

with 

ci^\t-T,Ax) = / dr'n'-^^t-r') - Ax. 



Above 0(2) denotes terms proportional to -^-^'P{x,t), m > 2. 

Denoting by L the characteristic length scale over which the probability 
distribution V{x, t) varies, we may introduce the new variable i = ^. Then 

rn^'t) = —-—Vii,t) 



Due to the dispersion of the probability distribution V{x,t) the character- 
istic length L asymptotically goes to infinity. 

Equating the coefficients of V{x,t) and -^Vixjt) on both sides of eq. 
lf2ijl . or equivalently, considering the 0th and 1st order in L, we end up 
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with 



dTK'^'^\t-T)z{T,t-T) = 1 (25) 
1 f°° 

- / dTK^^\t-T)z{T,t~T) = 



2 







dTK^i) {t-T)[[ dr'n^^'^ {t - t') -^\]z{T,t-T). (26) 
10 Jo 2 

Eq. lf26jl can be further simphfied using eq. lf25|l . which leads to 

) 

dTK^'^\t-T)z{T,t-T) = 

2 dTHl^^\t-T) dT'n^^\t-T')z{T,t-T)-l.{27) 



JO 



These two expressions, relating the asymptotic drift and diffusion proper- 
ties of a periodically driven renewal process as expressed by n'^^^t) and 
n'-^^t) to its microscopic properties defined by z{T,t) present the corner- 
stone result of our paper. Equations which govern the higher cumulant 
growth coefficients ^("-'(i), n > 3 can also be derived using this method by 
evaluating the coefficients of higher order derivatives of Vix, t) (for K,^'^\t) 

see appendix Ell- 
in a former work [5T| we investigated stochastic synchronization in a 
non Markovian two state model for excitable systems. The result presented 
in this paper for the mean velocity and effective diffusion coefficient can be 
shown to be equivalent to the results presented here (see appendix [Ejl . 

Often one is not directly interested in the number of events but in 
some quantity proportional to the number of events, Hke a phase, which 
increases by £ = 27r for each event, or a position which increases by some 
length £ = Z if we consider a unidirectional random walk with step-length 
I whose steps are governed by the described renewal dynamics. In these 
cases the ?th cumulant growth coefficient K^^\t) is scaled by £\ For K^^\t) 
and n'-^^t) this is achieved multiplying the constant term 1 on the right 
hand side of eqs. Ij25|l . and 112 7|l by £ or £^ respectively. 

Finally one may ask, why it is justified to prescribe a continuous Marko- 
vian envelope dynamics to an inherently non Markovian discrete process. 
The obvious idea, that the non Markovian nature of the discrete process is 
rendered Markovian by being mapped onto an extended continuous state 
space is misleading. The point is, that the continuous Markovian process 
x{t) as described by eq. Ijl6|l is not an envelope dynamics of the full dis- 
crete non Markovian process, but it only covers the asymptotic behavior of 
the non Markovian process, i.e the regime, where the discrete process has 
forgotten its initial preparation. 
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6 Comparison with known results for undriven 
renewal processes and periodically driven 
rate processes 

Let us evaluate expressions ll25|l and lf26jl for an undriven renewal process, 
i.e. z{T,t) = z{t). Then it follows that ^(^^(t) is constant and eq. Ij25|l 
leads to 

with 

poo 

(T") / dTr"w(r). 
Jo 

To derive this result we have used the fact that 

> poo n+1 1 

drr-zir) = / dr-— u;(r) = — - (T"+1) 
Jn n + 1 n + 1 ' 



which holds if z(r) decreases sufficiently fast for r ^ oo. Accordingly eq. 
l(27jl gives 

{Tf 

which agrees with the known results for stationary renewal processes 
The corresponding expression for K^^^{t) is presented in appendix IdI 
Next we consider a periodically driven Markov process, i.e. 

l-t + T 

w{t, t) = 7(i + t) exp ( — / dr'^^r')) 

and 

rt+T 

z{T,t) = c^p{^ dT'7(r')). 

Then it can be easily shown that eq. lf25|l is solved by 

^W(t)=£7W. 

The first term on the right hand side og eq. II27II can also be simplified in 
this case using integration by parts to give Therefore n^^^ (t) is governed 

by 



poo 

/ dTK^^'>{t-T)z{T,t-T) = £2 

Jo 
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which is solved by 

For more complicated processes with general time dependent waiting 
time distributions eqs. II25|I and 112 711 can only be solved numerically for 
the periodic solution. 



7 Numerical solution in Fourier space 

As eqs. Ij25|l and Ij26|l cannot be solved analytically for arbitrary wait- 
ing time distributions w(t, t — r) and corresponding survival probabilities 
z{t, t — r) one has to resort to numerical methods. To this end we perform 
a Fourier expansion of the periodic function n'^^^t) and analogously for 

oo T 



where Vl — 27r/T is the frequency of the external driving. We further 
expand the survival probability z{T,t) with respect to the second periodic 
argument as 

oo I 

z{T,t) — Zfe(r) exp(ifcf2t) Zkir) — — / dtz{T,t) exTp{—ikilt) 

Abbreviating 

/•oo /-oo 

Zk,i = / dTZfe(T)exp(-iZr2r) , hk,i ^ / dTTZk{T)exp{-ilV,T) {29) 

Jo JQ 

eq. Ij25ll can be written as 

C30 

K.'l!^^ Zm-k,Tn ^ S,S„i,0, m = -OO, . . . , CX) (30) 

k— — oo 

while eq. <t27ll reads 

oo oo 

^ ^ Zm — k.rn — 2 ^ ^ (^1) 

/c^ — oo /c^ — oo 

oo (1) (1) 

r '^z ^fc / \] (1) (1) J 
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The corresponding equation for the third coefficient k^^^ (t) in Fourier space 
is presented in appendix IdI 

These infinite dimensional inhomogeneous Hnear equations can then be 
numerically solved, after being truncated to a finite dimensional system. 

8 A simple example-comparison between the- 
ory and simulations 

Consider the toy renewal process, where the time between subsequent 
events is composed of a fixed waiting time, which depends on the signal 
phase of the previous event and a rate process with rate 7. The waiting 
time distribution is thus given by 

w{T,t) = 9{t - T(i))7e-'^(^-'^(*». 

Suppose further that the fixed waiting time is either Tq or Ti depending on 
whether the signal phase of the previous event was within [0, tt) or [tt, 27r), 
i..e 

^ ' \ Ti if m mod 2tt £ [tt, 27r) 

(32) 

A sketch of this system is shown in Fig. Elwhile corresponding waiting time 
distribution is plotted in Fig. 13 

Ti 




Figure 2: Depending on whether the periodic signal is in the first or second 
half period the system either waits the fixed time Tq or Ti . In both cases 
the system waits an additional exponentially with rate 7 distributed time. 

The corresponding Fourier coefficients Zk^i and hk,i as presented in eqs. 
H29|l can be analytically evaluated for this waiting time distribution, how- 
ever the final results, being too long and at the same time yielding not much 
information, will not be presented here. Having evaluated these Fourier 
coefficients, we calculated the mean frequency and the effective diffusion 
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Figure 3: Waiting time distribution of the toy model eq. l(32|l with di- 
chotomic fixed waiting time from eq. Ij32|l for Tq = 0.5, Ti = 1.5 and 
7 = 2. Depending on the signal phase of the event, the system waits ei- 
ther a long (dashed line) or a short time (solid line) plus an exponentially 
distributed time until the next event. 

coefficient according to eqs. l|3(11l and ll3Tll using LAPACK to solve these 
linear equations. The results are compared to simulations of the renewal 
process in Figs01and|3 showing perfect agreement. 

The mean velocity k'-'^' — 'T So dtfi^-^'>{t) and effective diffusion coef- 
ficient = dtK^^^t) can be used to characterize stochastic syn- 
chronization of the process to the periodic signal. These synchronization 
regions are defined by a rational relation between system frequency and 
signal frequency and a minimum of the effective diffusion coefficient. Al- 
though this model shows minima of the effective diffusion coefficient as a 
function of the driving frequency Fig. 0) these minima do not correspond 
to frequency synchronization as the system frequency and the signal fre- 
quency are not proportional. This stands in contrast to a similar model 
with a fixed waiting time and a dichotomically periodically modulated rate, 
used to describe periodically driven excitable systems j2l]. In this model 
we found several different n : m synchronization regions. Also the full pe- 
riodically time dependent coefficients k'^^ (t) and k^'^^ (t) as determined by 
our theory (I25|l and 112 711 agree with results taken from simulations of the 
underlying renewal process Fig. El Interestingly the effective diffusion co- 
efficient becomes negative for some values of the signal phase. However this 
does not imply that the periodic driving can be used to concentrate an en- 
semble of these systems as the period averaged effective diffusion coefficient 
K^^-* is always positive. 
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2 4 6 8 10 



Figure 4: Comparison of the period averaged mean frequency n^^^ = 
Y Jq dtn^^'' (t) and the period averaged effective diffusion coefficients k*^^^ = 
y Jp^ dtn'^'^^ (t) where T = 27r/f2 is the period of the signal, for the toy model 
eq. l(32ll with To = 0.5, Ti — 1.5 and 7 = 2. The solid lines are results 
of the theory eqs. (l25|l and lf27jl . numerically evaluated according to eqs. 

and H^illl truncated to 40 coefficients, while the symbols are obtained 
from direct simulations of the driven renewal process. The straight lines 
in the upper plot indicate n : m relations between system frequency and 
signal frequency, i.e. frequency locking. Clearly the system does not show 
this behavior. 
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t 

Figure 5: Comparison of the mean frequency n'-^^t) and the effective 
diffusion coefficients K^'^\t) for the toy model eq. (l32ll with Tq = 0.5, 
Ti = 1.5 and 7 = 2 and fl — 1.7. The soUd lines are results of the theory 
eqs. (j25ll and Ij27|l . numerically evaluated according to eqs. ^IM)^ and lj^^l|l 
truncated to 20 coefficients, while the symbols are obtained from direct 
simulations of the driven renewal process. 
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9 Conclusion 



We have considered the drift and diffusion behavior of periodic renewal pro- 
cesses and presented a general theory to express the mean frequency and 
effective diffusion constant and also higher order cumulant growth coeffi- 
cients in terms of the underlying time dependent waiting time distribution. 
The results of this approach where analytically shown to coincide with 
known results in the undriven case. For the periodically driven situation 
we confirmed the results of the theory by numerical investigations. We 
also showed agreement with a different approach for a more restricted class 
of renewal processes presented in [21]. The mean frequency and effective 
diffusion coefficient may be used to quantify stochastic synchronization. 
Taking into account the large amount of systems whose dynamics can be 
modeled as renewal processes, we anticipate a widespread applicability of 
our results, ranging from synchronization in neurons to the investigation 
of transport properties in molecular motors. 



A Expansion of the probability density gov- 
erned by a Kramers- Moyal equation 

Our aim is to express the phase distribution 'P{x — Ax,t — t) in terms of 
V{x,t) and its derivatives with respect to x, d" /dx"'V{x,t). To this end 
we start by expanding V{x — Ax, t — r)) in a Taylor series around x and t, 



V{x- Ax,t-T) = > > , , — 7^V{x,t) 



n—0 m— 



To process the time derivatives we use the Kramers-Moyal equation Ijl6|l 
taking care of the explicit time dependence of K^^^{t) which leads to 



V{x - Ax,t-T) = V{x,t) 



E 

771 — 1 



d_ 

dx 



V{x,t) 



m— 1 
m — 1 



m— 1 



— fm-l\ a™-l-'K(l)(t) 9'-lK(l)(t) 



1=1 



I 



2-1-1 



dx"^ 



V{x,t) + 0{i). 



where 0(3) denotes third or higher derivatives of V{x,t) with respect to 
X. The sums containing the coefficients ^("^(t) in a linear way be further 
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evaluated, leading to 

m=l m=0 "'U 

Jo 



'0 

The last term can be simplified to give 

m— 1 



^ ml ^ \ I J at™-!-' at'-i 

^=2 i=l ^ ' 

^ (m + l + 2)! V / + 1 / 

;=o m=0 V ' ' / \ / 



Jo ^' ^^"^ ^ tit 



. .... -- ^! at' 

m=0 Z=0 



Jo Jo 



Thus we eventually arrive at 



'P{x-Ax,t-T) = V{x,t) + c[^\t-T,Ax)—P{x,t) 



+cl^^(t-r,Ax)— P(x,i) + 0(3) (33) 



where 



c^\t-T,Ax) = [ dr' K^'^\t - t') - Ax 
Jo 



(34) 



and 



(ii^\t-T,Ax) = ^-AxTdr'n^^^t-r') (: 
2 Jo 

-1 r dr'K^^\t-T')+ r dT'K^^\t-T') r dT"K^^\t-T") 

2 Jo Jo JO 
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B Relation between moments and cumulants 
of the number of events of a renewal process 
and its characteristic functional 



Consider the moments M^"^ defined by (c.f. eq ^) 



exp[^G.(e^-l)^] 



z=0 



(36) 



Generally the relation between moments and cumulants is given by 

5:i,MW.exp[^i,i^W]. (37) 



k=l 



From the relation 



fc=0 



fc! 



2=0 



and eq. l(3fijl we deduce (taking into account the analyticity at z = of the 
functions involved) 



C30 t. 



^ijMW=exp[5:G.(e^-l)^] 



s=l 



and thus according to eq. Ij37|l 



oo u 



EiT^''=^=EG.(e^-l)^ 



fc=i 



With 



dz" ^ fc! 



z=0 



the final result 



on 



follows. 
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C Relation between the Kramers- Moyal coef- 
ficient and the growth of the cumulants 

Consider the stochastic process x{t) governed by 

r\ OO / 1 \ 71 Pin, 

a^^(-'^) = E^«^"HO^^(-'^) (38) 

We are interested in the grows of the cumulants K^'"\t) of x{t). The 
moments 

/oo 
dxx^-Vix, t) 
-co 

obey 

Assuming further that V{x, t) decreases sufficiently fast for x —>■ ±oo such 
that 

lim x"-V{x, t)=0 

X— '■itoo 

the above expression can be evaluated using integration by parts to give 



Now the moments and cumulants are related by 

oo ^ oo ^ 

^|^MW(t)=exp[5:i^i^W(0]. (40) 
and thus by differentiating this equation with respect to t 

ooj^7 k ^ k 

k\ dt ^ ^ f-^ k\dt ^ ' f-' k\ ^ ' 

k=0 k=l fe=0 
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Inserting the moments dynamic eq. into the left hand side of this 

equation, it can be easily checked that 

D The third cumulant growth coefficient K^^\t) 

The third coefficient k*^'^^ (t) can be evaluated by equating the coefficients 
corresponding to the order ■^'P{x,t) in eq. (I24II . Inserting the expansion 
into this equation leads to 



1 

6 7o 



dTK(3)(i _ r)z(r,i - r) = (41) 

1 PCX. PT 1 

- / dr^(2) ^t-r)[ dr'n^^) [t - r') - -] z(r, < - r) 



+ / dT'K^^\t-T') / rfT"K(l)(<-T")]z(T,t-r) (42) 

Jo Jo 
or using eqs. lf25|l and ll27ll 

/"OC /"OO 

/ dTK^^\t-T)z{T,t-T) = / dTz{T,t-T) (43) 



3K(')(i-T) / dr'K(i)(t-r') + 3KW(i-^) / dr'^^^) (t - r') 
Jo Jo 

-6kW(*-^) C dT'K^^Ht-T') f dT"K^^\t~T") 



1 

For an undriven renewal process eq. (j43ll can be directly solved, leading to 
the time independent coefficient 



For a rate process with periodically modulated rate ^{t) l(i3| eventually 
leads to 

K(3)(t) = £3^(i). 
Generally, the periodic solution 

K(3)(i)= 4^^exp(ifcm), ^ = -r / dtK(3)(t)exp(-ifcrJt), 
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of eq. H'M can be numerically obtained in Fourier space as a solution of 
the (infinite) set of linear equations 



,(3) 

fc— — oo 



fc— — OO l — — OQ,l^0 



oo oo oo 

E [[ E E 



'6 

A.;^ — CXj l — — OQ.l^0 j — — OO,j^0, — l 

' (^~j^{^7n — k — l—j,m—j ^m — k—l—j,m—j — l) 



^fe '^Z '^J /I 



^" . ^ ^ {^^m—k—l—j.jn—j — l — fc — /—j.m )) ] 



oo (1) (1) (1) 

ic fcr" (i* 

"r / ^ /-^^^"^ \^va — k^m ^m — k,m—l < '^^^^'^rn — k^mj 
°° ^(1)^(1)^(1) 

V- ^fc ^0 (u , , , , ^. ^(1)^(1)^(1), , 

where we used the Fourier decompositions of kM'' and k^"^^ according to eqs. 

Ij28|l while z^^i and hk,i are defined in eq. (I29II and jk,i is defined as 



Jk,l 



1 f°° 

2 Jo 



For the toy model introduced in section |S| we have evaluated the period 
average of n^^^t) both according to the theory, eqs. and ijiiji and 

from direct simulations of the driven renewal process. Both results agree 
very well (Fig. ©), thus confirming our theory. 

E Equivalence with the two state model for 
excitable systems 

In a former work [21J we investigated stochastic synchronization in a non 
Markovian two state model for excitable systems. Such an excitable dy- 
namics can be approximated by two steps, namely an excitation from the 
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Figure 6: Comparison of k^^^ — t Io dthS^\t), where T = 2it/Q, is the 
period of the signal, for the toy model eq. Ij32|l with Tq = 0.5, Ti = 1.5 
and 7 — 1.7. The solid lines are results of the theory eq. numerically 
evaluated according to eq. ijiiji truncated to 20 coefficients, while the 
symbols are obtained from direct simulations of the driven renewal process. 
Both results agree very well. 

stable fixed point (rest state), which can be well modeled by a rate process, 
followed by a long excursion in phase space, which constitutes the excited 
and refractory state. This excursion takes some time which is distributed 
according to a waiting time distribution. As the excitation step is much 
more sensible to the periodic driving than the time spent on the excitation 
loop, we may assume this waiting time distribution w{t) to be independent 
on the running time. The rate 7(t) of the excitation step however, depends 
on the periodic signal and thus periodically on the time t. 

In this model we have defined a phase in order to investigate synchro- 
nization between the signal and the system's dynamics. This phase was 
chosen to increase by 27r each time the system has been excited and re- 
turned back to the rest state. We have derived expressions for the mean 
phase velocity u){t) and effective phase diffusion constant D{t). In the fol- 
lowing we show that these expressions agree with the corresponding expres- 
sions of our general theory presented here. Namely by choosing £ = 27r the 
phase velocity can be identified with K^^\t) while the effective phase 
diffusion coefficient D{t) in this prior work differed from our definition of 
K^^\t) by a factor of 2, i.e. n'^'^\t) = 2D{t). 

In this two state model the mean phase velocity was given by 

c^W = 2^7(i)<zf (i) (45) 
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where the dynamics of ^2^'' (t) are governed by 

1 - it) = / dr^it - T)gf {t ~ t)S{t) (46) 
Jo 

The effective diffusion coefficient was calculated as 

Dit) = -27rj{t)qi^\t)+7roj{t) (47) 



with 



POO 

2T,q^'> it) = 27r / dTS{T)jit - T)qi^'> [t - r) (48) 
Jo 

/•OO /"T 

+ / dTz{T)uj{t ~ t){ I dT'uj{t-T')-2TT). 



The auxiliary variables q^^ and 92^'' can be eliminated leading to 

/■OO 

7(0 / ci™(t - T)z(r) + u}{t) = 27r7(i) (49) 
Jo 

and 

/•OO 

7(i) / dTD(t-T)z(T) + D(t) = (50) 
Jo 

/"OO /"T 

-27r^7(t) + 2T:w{t) + '-f{t) / dTUj{t ~ t) druiit - t')z(t). 



In the following we show that these equations for w(i) and D{t) are equiv- 
alent to our general eqs. Ij25|l and Ij27ll which for w(i) and D{t) read 



dTUj{t-T)z{T,t-T) = 2ti (51) 

(52) 

and 



nOO 

I dTD{t ~ T)z{T,t - t) = 

Jo 



dTUj(t-T) / dT'w(i - T + r')z(r, i - r) - 27r2 



(53) 



The waiting time between two events is given by the sum of the excitation 
time and the time spent on the excitation loop. The time dependent waiting 
time distribution for one cycle is therefore given by 

t+r 

'I .J—'l\ 



w{T,t) = I dT'w{T')-f{t + T)exp{~ / dT"7(T")). 

Jt + T' 



(54) 
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while the corresponding survival probability reads 

z{T,t) = / dT'w(T') exp(— / dT"'^{T")) + z{t) (55) 

Jo Jt + T' 

where z{t) = 1 — /J" dT'-wlr) is the probability to spent a time longer than 
r on the excitation loop. From these equations we obtain the relation 

7(t)z(T,i-r) = iy(T,t-r) +z(t). (56) 

Multiplying eqs. ll5T|l and l|53ll by 7(i) and replacing the time dependent 
survival probability z(t, t — r) according to eq. ll56|l we obtain 



7(i) / dTLu{t ~t)S{t) + dTuj{t-T)w{T,t~T) = 27r7(t) (57) 
Jo Jo 

and 

7(i) / dT[D{t - t) - Lj{t ~ t) / dT'uj{t-T-T')]S{T)+ (58) 
JO Jo 

dr [D{t ~ t) - uj{t - t) [ dT'uj{t - t - t')] w{t, t - t) = -2Tr^j{t) 



The second term on the left hand side in both equations can be further 
simplified. To this end we differentiate eqs. ll5Tll and lf53jl with respect to t 
(see appendix EJ; using eq. 115 111 to simplify the time derivative of eq. Ij5^i|l 
which results in 



uj{t) ^ / dTUj{t - t)w{t, t - t) (59) 
Jo 

and 

/"OO pT 

D{t) ^ / dT[D{t - t) - L0[t - t) dT'u;{t-T-T')]'w{T,t-T)+2TTUj{t). 

Jo Jo 

(60) 

Note that the above equations can be equivalently derived from eq. lf23jl 
by performing the same procedure we have applied to eq. Ij22|l in order to 
obtain eqs. (ESJ and lf?7jl . 

Inserting finally eqs. Il59|l and Ij60ll into eqs. Ij57ll and Ij58|l . eventually 
leads to the two state model eqs. ijiflll and l|5(11l . 

Thus we have shown that solutions of D{t) and uj{t) of the general 
theory are also solutions of the two state model for excitable systems pre- 
sented in (21]. As these solutions are in general unique we have shown the 
equivalence between the two models. In pi] these solutions were compared 
with simulations of the renewal process, showing total agreement within 
simulation precision. 
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F Time derivatives of integrals involving time 
dependent survival probabilities 

Let w(r, t) be a time dependent waiting time distribution, i.e. 



w{t, z) >Q V T,t and 
Let further 



dTw{T, t) = 1 V r. 



z{T,t) = 1 - / dT'w{T',t) 

Jo 

be the corresponding survival probabiUty. Then for any sufSciently well 
behaved function g{t) we have 

d f°° 

- drgit - t)z{t, t-T)= g{t) - dTg{t - t)w{t, t - r). (61) 

This can be seen as follows: 
d f°° 

di Jo '^^^^^ ~ ~ 
= j dT(^-^g{t-T)z{T,t-T)+g{t-T)^z{T',t-T)l,^^) 

= -g{t - r)z{T, t-r) - dTg{t - t)w{t, t - r) 
-r=o Jo 

POO 

= 9{t) - / drgit - t)w{t, t - r). 
Jo 

where in the last step we assumed that g{t — t)z{t, t — r) converges to zero 
as T — »■ 00. 

We further have the identity 

d f°° f 

- / dTg{t -t) I dr'git - r + t')z{t, t - r) 

= dri^-— [g{t - r) dr'git - r + t')z{t, t - r) 

+g{t - T)g{t)z{T, t-T)-g{t-T) dr'git - r + r')w{r, t - r)) 

POO 

= g{t) / drg{t - r)z{r,t - r) (62) 
Jo 

poo pT 

- I drg{t-r) / dr' g{t - r + r')w{r,t - r). 
Jo Jo 
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